﻿#pragma once

#include <xy/solver/linear.h>
#include <xy/solver/nonlinear.h>

namespace xy
{

namespace solver
{

using time_grid = std::vector<math::vecxd>; // 迭代网格，按照时间对应向量
using gtimer = time_grid::iterator;         // 网格迭代器

struct integrator
{
    std::size_t series; // 级数
    double currentTime; // 当前时间
    double stepTime;    // 时间步长

    integrator(std::size_t series, double stepTime);

    virtual gtimer shift(funcxxt, gtimer); // 时间算子作用，返回作用后的时间
};

// 多步法
struct multi_step : integrator
{
    math::vecxd alpha; // u 系数，s s-1 s-2 ...
    math::vecxd beta;  // f 系数，s s-1 s-2 ...

    multi_step(std::size_t series, double stepTime);

    /**
     * @brief 通用隐式多步法
     *
     * @param[in] f
     * @param[in] U
     * @return gtimer
     */
    gtimer shift(funcxxt f, gtimer U) override;
};

struct adams_bashforth : multi_step
{
    adams_bashforth(std::size_t series, double stepTime);

    /**
     * @brief 显式多步法
     *
     * @param[in] f
     * @param[in] U
     * @return gtimer
     */
    gtimer shift(funcxxt f, gtimer U) override;
};

struct adams_moulton : multi_step
{
    adams_moulton(std::size_t series, double stepTime);
};

struct bdfs : multi_step
{
    bdfs(std::size_t series, double stepTime);
};

// 单步法
struct one_step : integrator
{
    math::matxd A;       // 系数矩阵
    math::vecxd weights; // 加权
    math::vecxd nodes;   // 节点
    const std::size_t p; // 阶数，用于外插法

    one_step(std::size_t series, double stepTime, std::size_t p);

    /**
     * @brief 理查德外插法，p 阶方法精确到 p+1 阶
     *
     * @param[in] f
     * @param[in] U
     * @return gtimer
     */
    gtimer extrapolate(funcxxt f, gtimer U);
};

struct ERK : one_step
{
    ERK(std::size_t series, double stepTime, std::size_t p);

    /**
     * @brief 时间迭代
     *
     * @param[in] f
     * @param[in] U
     * @return gtimer
     */
    gtimer shift(funcxxt f, gtimer U) override;

    /**
     * @brief 迭代到下一步之前，计算出方程组的根
     *
     * @param[in] f
     * @param[in] U
     * @return std::vector<math::vecxd>
     */
    std::vector<math::vecxd> preshift(funcxxt f, gtimer U);
};

struct IRK : one_step
{
    IRK(std::size_t series, double stepTime, std::size_t p);

    /**
     * @brief 时间迭代
     *
     * @param[in] f
     * @param[in] U
     * @return gtimer
     */
    gtimer shift(funcxxt f, gtimer U) override;

    /**
     * @brief 迭代到下一步之前，计算出方程组的根
     *
     * @param[in] f
     * @param[in] U
     * @return math::vecxd
     */
    math::vecxd preshift(funcxxt f, gtimer U);
};

struct classical_RK4 : ERK
{
    classical_RK4(double k);
};

struct ESDIRK6 : IRK
{
    /**
     * @file RungeKutta.h
     * @author xingyifan
     * @date 2023-04-24 20:28
     *
     * @description: 通过时间步长构造
     */

    ESDIRK6(double k);

    /**
     * @file RungeKutta.h
     * @author xingyifan
     * @date 2022-07-28 22:16
     *
     * @description: 对角隐式方程，当函数为线性形式：f(X,t) = BX + g(t), B = I - k a_ii A 时，可以通过求解线性系统推进
     * 不能单独求解每个方程，因为它们互相关联，需要同时求解
     */

    gtimer shift(const math::matxd &B, func1x g, linear_ptlu &solver, gtimer U);
};

struct gauss_legendre_RK : IRK
{
    /**
     * @file RungeKutta.h
     * @author xingyifan
     * @date 2023-04-24 20:28
     *
     * @description: 通过阶数和时间步长构造
     */

    gauss_legendre_RK(std::size_t s, double k);

    /**
     * @file RungeKutta.h
     * @author xingyifan
     * @date 2022-07-27 23:01
     *
     * @description: 纯隐式方法，使用向量牛顿迭代法求解
     */

    gtimer shift(funcxxt f, gtimer U) override;
};

struct embedded_RK : ERK
{
    math::vecxd dualWeights; // 对偶 RK 权
    const std::size_t dp;    // 对偶精度

    embedded_RK(std::size_t series, double stepTime, std::size_t p, std::size_t dp);

    /**
     * @brief 自适应步长（暂时不测试）
     *
     * @param[in] f
     * @param[in] U
     * @return gtimer
     */
    gtimer auto_step(funcxxt f, gtimer U);
};

struct fehlberg6 : embedded_RK
{
    fehlberg6(double k);
};

struct dormand_prince7 : embedded_RK
{
    dormand_prince7(double k);
};

} // namespace solver

} // namespace xy
